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ABSTRACT 

Radiative feedback is among the most important consequences of clustered star for- 
mation inside molecular clouds. At the onset of star formation, radiation from mas- 
sive stars heats the surrounding gas, which suppresses the formation of many low- 
mass stars. When simulating pre-main-sequence stars, their stellar properties must 
be defined by a prestellar model. Different approaches to prestellar modeling may 
yield quantitatively different results. In this paper, we compare two existing prestellar 
models under identical initial conditions to gauge whether the choice of model has 
any significant effects on the final population of stars. The first model treats stellar 
radii and luminosities with a ZAMS model, while separately estimating the accretion 
luminosity by interpolating to published prestellar tracks. The second, more accurate 
prestellar model self-consistently evolves the radius and luminosity of each star under 
highly variable accretion conditions. Each is coupled to a ray tracing-based radiative 
feedback code that also treats ionization. The impact of the self-consistent model is 
less ionizing radiation and less heating during the early stages of star formation. This 
may affect final mass distributions. We noted a peak stellar mass reduced by 8% from 
47.3M0 to 43.5M0 in the evolutionary model, relative to the track-fit model. Also, the 
difference in mass between the two largest stars in each case is reduced from IAMq to 
7.5Mq. The HII regions produced by these massive stars were also seen to flicker on 
timescales down to the limit imposed by our timestep (< 560 years), rapidly changing 
in size and shape, confirming previous cluster simulations using ZAMS-based estimates 
for prestellar ionizing flux. 

Key words: hydrodynamics - radiative transfer - stars: formation - stars:pre-main- 
sequence - stariprotostars - ISM: HII regions. 



1 INTRODUCTION 

The conversion of molecular gas into fully-formed stars is 
complex, involving several diverse processes. These different 
processes are linked to each other through feedback mecha- 
nisms that make isolating and understanding the contribu- 
tion of each process a difficult task. A key point in this re- 
gards is that stars also rarely form in isolation, but instead 
are seen to be fo rming in clusters and subclusters wi thin 
molecular clouds (jClarke et al.ll200Ql : ll^sti et al.ll2000l ). In 
the cluster environment, the formation of a sufficiently mas- 
sive star can affect all the others through the energy it ra- 
diates back into the cloud. Numerical simulations of star 
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formation have made it very clear that the effects of stellar 
radiation cannot be neglected. Simulations including some 
form of radiative transfer show a dramatic reduction in 
the production of brown dwarfs and other low-mass stars 
(Offner et al. 2009), due to an increase in ^as t emperatures 
reducing fragmentation (|Krumholz et al.ll2007l : IPeters et al.l 
l2010bl ). More of the available gas mass ends up being ac- 
creted by the fewer, larger stars formed, and the fragmen- 
tation that do es occur takes place i n optically thick s elf- 
shielding discs (jKrumholz et al.|[2007l : IPeters et al ] |2010al lbh. 
The fact that radiation affects the mass spectrum in simu- 
lations of molecular cloud clumps has obvious implications 
for the shape of the initial mass function, for example the 
suppressio n of excessive brown dwarf form ation ( B atell2009l : 
iKrumholz et al.ll2010l : [Peters et al.ll2010bl ). 



© 2011 RAS 



2 M. Klassen, R.E. Pudntz, & T. Peters 



Massive stars al so emit prodigious amoun ts of UV radi- 
ation (Hoa re et al.. „2007l : iBeuther et alj |2007l) creating ex- 
panding HII regions. The hot (10^^ K) gas expands into 
the colder (10^ K) surrounding low-pressure gas, creating 
another feedback mechanism and ionized region that may 
contr i bute t o the destruction of molecular clouds ( | Ketol 
l200d, l2003l . l2007l : iMatznerl l2QQ2l : IPeters et alJ l201Q allbh. 
HII regio ns can be observed by th eir radio continuum 
emission (|Mezger Henderson 1 19671), or by the ir recom- 
bination lines (e.g. Wood & Churchwelll (|l989l ) use the 
H76a line). More recently, observ ations have shown time 
varia b ility in HII regions ([Franco- Hernand ez Rodriguez. 



20041: iRodrfguez et al.1 l2QQ7l: iGal van-Mad rid et al.1 



Gomez et al.l l2008l V iFranco- He rnandez & Rodriguez 
have suggested that such observed time-variability may be 
due to the changes occurring in the source of the ionizing ra- 
diation, though it may also be due to in creased absorption in 
the rapidly-evolving core of the nebula. IPeters et al. 1 (l2010d l 
present a technique for using synthetic radio maps to study 
the time-evolution of stars forming in a cluster environment 
and variability in the morphology and size of HII reg ions. 
Analysis of these simulatio ns by iPeters et al.l (|2010al l and 
iGalvan-Madrid et al.l (|201lh confirmed variability in the flux 
and size measurements of HII regions, which in a few cases 
might be observable on timescales of ^ 10 years. They also 
noted that positive changes were more likely to occur than 
negative changes, i.e. that most of the flux variations were 
increases rather than decreases. 

To further explore the impact of radiative feed- 
back and the possible variability in HII regions, simula- 
tions must be equipped with go od protostellar r aodels . 
These have been investigated by Palla & Stahler ('iQQlJ), 
Stahler (1992), Nakano et al. (2000), McKee & T^ 
, Offner et al.l ([20091) and Hosokawa & Omuk3 
, among others. It is clear from these models, that 
the evolution of a protostar depends heavily on the mass ac- 
cretion rate. Among other things, they show that the radius 
of the protostar may grow or contract depending on the 
stellar evolutionary stage. With a radius that can change 
significantly during the pre-main-sequence lifetime of the 
star, the effective temperature can also be expected to vary 
significantly. To study this, we simulate the formation of 
a cluster of stars inside a molecular cloud. We equip the 
stars with o ne of two prestellar mode ls based on the on es 
described in IPeters et aD (|2010al l and I Offner et al.1 (120091) . 
each with its own characteristics. The Offner et al.l (12009 ) 
model has already been used to study star cluster forma- 
tion in Krumholz et all (|201lh . though with different initial 
conditions. Ours is the first simulation with the protostellar 
model to also include the effects of ionizing radiation and 
HII region formation. We connect the model to a radiative 
transfer method that computes the heating and ionization 
due to radiation from the stars formed in the simulation. 

The differences between the two m odels is expla i ned in 
[2Tl but the key difference is that the lOffner et all (|2009l ) 
model treats the evolution of the radius and luminosity self- 
consistently. The choice of stellar model affects the early 
evolution of stars in a cluster, and may have repercusions 
for the final mass spectrum. Though not entirely conclu- 
sive, we find that reduced heating and ioniza tion in the 
early stages of star formation when using the Offner et al.l 
(|2009l l model resulted in a more equitable mass distribution. 



With the"PeterseiaL| (|2010a') model, the cluster came to be 
dominated more by a single star about IAMq more massive 
than the nex t largest, compared to a 7.5M0 gap in the 
lOffner et al.l ([200s) model simulations. 

Other effects of the self-consistent prestellar modeling 
are delayed ionization of the cluster gas by 3% of a freefall 
time (17.7 kyr), and delayed heating of the cluster gas by 
1% of a freefall time (5.9 kyr). 

Our numerical approach is described in Section [2l In 
Section [3] we list our results for the early evolution of star 
clusters with massive stars. Our assessment of the impact of 
protostellar modeling we discuss in Section[3]and summarize 
our findings in Section [5] with a view to future simulations. 



2 NUMERICAL METHODS 

We perform nu merical simulations using the FLASH hydro- 
dynamics code (|Frvxell et al .11200^ 1 in its version 2.5. It is an 
adaptive-mesh refinement code that solves the gas-dynamic 
equations on an Eulerian grid and inclu des self- gravity, 
cooling by dust and by molecular lines (Baneriee e t al I 
I2OO6I ). and radiative transfer. It has been modified to 
include Lagrangian sink particles (|Baneriee et al.l l2009l : 
iFederrath et al.ll2010l l to represent (proto) stars, and a ray- 
tracing scheme to handle ionizing and nonionizing radiation 



feedback from stars originally developed bv lRiikhorst et al 
(.2006i). then extended and optimized by IPeters et al 
(|2010al l. They also testing the code against handling a D- 
type ionization front, comparing it to the approximate so- 
lution found by Spitzer (1978), while the code's ability to 
handle R-t ype ionization fronts has already been tested by 
llliev et aP (|2006 ). Accretion rates onto sink particles are cal- 
culated based on a single time step. FLASH does not have 
adaptive time steps, so every refinement level advances with 
the same time step. 

The opacities for the no n-ionizing radiation are the 
same as in Peters et al. I (|2010a|V, We use Planck me an opaci- 
ties as interpola t ed fro m the IPollack et all (| 19941 ) data by 
iKrumholz etHI (|2007l ). They assume that the radiation 
temperature is equal to the gas temperature because their 
core is optically thick. We make the same approximation 
using the assumption that the star will be embedded in an 
(unresolved) dense envelope of gas through which the stellar 
radiation must propagate before entering the scales of our 
simulation, thereby changing its spectrum accordingly. 

We subsequently added an additional module to han- 
dle the protostellar evolution of our sink particles, which 
is based on a subgr id physics model described in detail 
m lOffner et al.l (|2009l V The protostellar model connects di- 
rectly to the radiation module so that stellar surface tem- 
peratures and stellar radii are handled self-consistently. The 
physical stages of this model are listed in Table [21 



2.1 Protostellar models 

The radiative feedback model is coupled directly to the sink 
particles. Rays are cast outwards from each sink particle 
and the column density along each ray computed usin g the 
hybri d-characteristics scheme described in Riikhorst e t al.l 
((20061). At each cell in the computational domain, the pho- 
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toionization rate and heating rate are calculated. These are 
set by the specific mean intensity along the ray, 
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which depends on a knowledge of the radius of the star rg^Q^^.. 

The stellar radius depends on the choice of stellar 
model. The simplest stellar model would be to assume that 
all stars are ZAMS stars and use a lookup table, such the 
one by Paxton (2004), to retrieve the radius for a star re- 
siding in a particular mass bin of the table. Such a table 
will also contain surface effective temperatures for ZAMS 
stars. The intrinsic stellar luminosity is then found from 
^int ~ 4:7TRlaT^^. This model may be acceptable for most 
circumstances, but breaks down for when attempting to 
model pre-main-sequence stars. If one treats these low-mass 
stars as ZAMS stars, the model will underestimate their 
radii and overestimate their surface temperatures. It will 
also lead to an overestimation of the accretion luminosity 
Lace = GMM/R. 

In IPeters et al.l (|2010al l. a kind of "augmented ZAMS" 
model is used, which we'll refer to as A-ZAMS throughout 
the paper. This prestellar model uses a ZAMS description 
as detailed above when calculating the stellar radius and 
instrinsic luminosity of stars. To avoid overestimating the 
accretion luminosity, a separate accretion radius is calcu- 
lated. This is achieve d by referencin g the pre-main-sequence 
tracks computed by iHosokawa Omukai (2009) for mass 
accretion rates between lO^^Mo/yr and 10~^ Mq/jy and 
then interpolating between them based on the current mass 
accretion rate for the star. The advantage of this model is 
that it is relatively straightforward to implement and pre- 
vents grossly overestimating the accretion luminosity, which 
dominates the total luminosity of a star during its early life- 
time. The disadvantage of this model is that it is not self- 
consistent and relies on two separate radii being computed 
or retrieved from a table. The accretion radius, found by in- 
terpolation to tracks of constant accretion rate, is sensitive 
to fluctuations the accretion rate. A rapidly fluctuating ac- 
cretion rate means the accretion radius will fluctuate with 
equal rapidity — and unphysical consequence. 

An alternative approach is to use the self- consistent 
evolvi ng protostellar model de veloped by ll^n McKee 
(HqqI) and described in detail in lOffner et aD (|2009l ). Stars 
are modeled as polytropes and every sink particle in our 
simulation is assigned several additional properties: a stel- 
lar radius rg^ar? an intrinsic luminosity L-^^^, a polytropic 
index n, an unburned deuterium mass m^, and a nuclear 
burning evolutionary stage. At every timestep in our simu- 
lation, we evolve this ha ndful o f variables according to the 
equations given in Offner et al. I dloog). The model is based 
on a one-zone protostellar evolution model introduced by 
iNakano et al.l ([l995i) and f u rther developed bv lNakano et al.l 
(|200Q| ) and lTan k McKe'3 (|2004l l. 

We refer to this prestellar model as the "evolving pro- 
tostar" model, to distinguish it from the A-ZAMS model 



employed in IPeters et al.l (|2010ah and used for comparison 
here. It is so called because the stellar properties are co- 
evolved with the rest of the simulation instead of calculated 
on-the-fly. 

When a sink particle's mass exceeds O.IM©, we activate 
our protostellar evolution code and initialize the radius and 
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Figure 1. Luminosity evolution for protostars accreting mass at 
various rates. The solid lines show the intrinsic (stellar) luminos- 
ity following the evolving protostar model, whereas the dashed 
lin es show the luminosity deriv ed from stellar structure modeling 
bv lHosokawa h OmuM (|2009l) . 



polytropic index respectively as 

Am/At 



2.5 K 







n = 5 — 3 



1.475 + 0.07 log 



10-5 Mo yr-i 
Am/At 



(2) 
(3) 



For cool stars, the Hayashi limit sets the luminosity, but 
above this a main sequence luminosity is assumed. Thus, 
^int ~ max(Lif , Lms), with Lh — 4:7TR^aT^ and Th = 
3000 K. ZAMS values for the radius and lumino s ity ar e com- 
puted using the fitting formulas bv lTout et aD (1 19961 ). 

Apart from initializing our model at a higher starting 
mass, the only other significant difference is that we take 
accretion luminosity to be Lace = GMM/R. The evolv- 
ing protostar model treats accretion onto the disc, with an 
associated luminosity of L^^g^ = (1/2) GMM/R (standard 
for an alpha disc), and surface accretion, with luminosity 
Lace = (1/4) GMM/R. This is due to the assumption that 
some of the energy is being used to drive a wind. We do not 
make this assumption. We also rely on tables of polytropic 
stellar parameters that we computed ourselves. In all other 
respects, our protostellar model follows the one described in 
lOffner et al.l ((2009|). 

Protostars evolve through multiple distinct nuclear 
stages in this code during which the radius is at times 
expanding (such as during the early accretion phase) and 
at times contracting (such as during the end stage as the 
protostar approaches the main sequence to be a mature 
star). Once our stars reach the main sequence, we assign 
them a radius and lum inosity based on the fitting formu- 
las of iTout et a l.' (1996). We neglect any special treatment 
of metallicity-related effects and consider only stars of solar 
metallicity. 
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Table 1. Runtime parameters of the clustered star formation simulations 



Run 


Mass 


Density profile 


Temp 


Rotation 


Stellar Model 


Feedback 




1 


1OOOM0 


^-3/2 


30 K 


13 = 0.05 


"Evolving Protostars" 


Radiative; 


raytracing method 


2 


1OOOM0 


^-3/2 


30 K 


^ = 0.05 


"Augmented ZAMS" 


Radiative; 


raytracing method 
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Figure 2. Mass-radius relation for accreting protostars. The 
solid lines show the stellar radius following the evolving proto- 
star model, whereas the dashed lines sh ow the radius derive d 
from stellar structure modeling by Hosoka wa &: Omuka 

1 (!2009lV 

The one-zone model approximates the stellar structure results to 
within a factor of ~ 2. 

To compute the ionizing flux, we take the stellar radius 
and surface temperature from either a table of ZAMS val- 
ues (in the case of the A-ZAMS model) or read the current, 
evolved values from the sink particle properties computed 
by the protostellar code (in the case of the evolving proto- 
star model). The flux of ionizing photons is the computed 
by integrating the Planck function above the threshold fre- 
quency for hydrogen ionization. The radiative feedback code 
computes gas heating considering both the intrinsic and ac- 
cretion luminosities. 

This protostellar evolution code is a one-zone model 
that upgrades the current treatment of sink particles in 
FLASH and is a more accurate representation of pre-main- 
sequence stars. In Figures [1] and [2] we compare the re- 
sults of this model with th e stellar structure modeling of 
iHosokawa Omukail (|2009l ). which is expected to be more 
accurate than one-zone modeling. 

We compare the behaviour of our code at different ac- 
cretion rates ranging between a slow lO^^Mo/yr to a rapid 
lO^^Mo/yr. These represent the typical range of accretion 
rates we see in our simulations and expect of stars form- 
ing in clusters within molecular clouds. The stability of the 
code was tested over a range of accretion rates and timestep 
sizes. Although our tracks do not agree perfectly with the 



iHosokawa Omukail (|2009l l simulations, the agreement is 
to within a factor of ~ 2. 

Table [21 with its accompanying figure, sum marizes what 
i s des cribed in detail in the appendices of lOffner et al. 
(|2009h . 



2.2 Initial conditions 

The strength of the radiative feedback code we employ lies 
in its ability to produce realistic HII regions. It was be- 
lieved, however, that since the radius and stellar luminosity 
of you n^ protostars are rep resented by ZAMS-equivalent val- 
ues m IPeters et all (|2010al l. that the ionizine: flux would be 
overestimated. To study whether this was indeed the case, 
and also what impact a different prestellar model would have 
generally, we simulated a collapsing molecular cloud clump 
with each of the two prestellar models. 

In the first case, we chose to repeat the cluster simula- 
tions described in IPeters et al.l ((2010a|) with similar initial 
conditions, but at a slightly lower resolution. Because we use 
the same FLASH code, sink particles, and radiative feedback 
code, we can isolate the effect of including a protostellar 
evolution model. We begin with a IOOOM0 self- gravitating 
clump of molecular hydrogen at an initial temperature of 
30 K. The cloud is in solid body rotation with a ratio of 
rotational to gravitational energy of /3 = 0.05. Our simula- 
tion box is 3.89 pc on a side. At maximal refinement, the 
grid size is 196 AU. The density profile features a flat cen- 
tral region extending out to a radius of 0.5 pc, then falling 
off according to an t~^^^ power law. The central density 
is pc = 1.27 X 10~^° g cm~^. The density drops off until 
reaching an ambient cutoff density of pg^t ~ ^-^^ ^ 10~^^ 
g cm""^. Sink particles have a radius of 1175 AU, or 6 times 
the grid size at maximal refinement. The cut-off density for 
sink particle creation is 4.4 x 10~^^ g cm~^. 

We note that clumps of t h is size and niass are expected 
to be turbulent (|Blitzf Il993l : lEvand Il999l : IWilliams etHI 
I2OOQI ). However, in order to build up physical understand- 
ing of the complex process of cluster formation, we follow 
IPeters et al.l (|2010al ) in this study and ignore turbulence so 
that we can isolate the important radiative feedback effects. 
Turbulence will be added in subsequent papers. 

We show the results of two of our simulations: one using 
the A-ZAMS approach for stellar effective temperature and 
stellar radius, and a second with the evolving protostellar 
model. 

As the simulation progresses, the original mass profile of 
the clump quickly disappears as it undergoes gravitational 
collapse to produce a rotating central disc. Stars, repre- 
sented by sink particles, are allowed to for m when the local 
condit ions satisfy the criteria described in iFederrath et al.l 
(|201Q| ). 
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Table 2. Description of the stellar evolutionary stages in the evolving protostar model 



Stage Features 






Pre-Collapse 




Mass m < 0.1 Mq 

Cannot dissociate H2 and cause second collapse to stellar densities. 


1 


No Burning 




Object has collapsed to stellar densities. 
Tc still too cold to burn D. 
Tc < 1.5 X 10^ K 

Radiation comes purely from gravitational contraction. 
Star is imperfectly convective. 


2 


Core D burning 


at fixed Tc 


Temperature reaches required Tc ~ 10^ K to burn deuterium. 
D burning acts as a thermostat keeping temperature constant. 
Star is fully convective. 


3 


Core D burning 


at variable Tc 


D is exhausted. 

Core temperature now rising again. 

Star remains fully convective. 

Accreted D dragged down to core and burned. 

Rising core temperature reduces opacity. 

Convection in the stellar core eventually shuts down. 


4 


Shell D burning 




Star core changes to a radiative structure, swelling the radius. 
D burns in a shell around the core. 

After initial swelling, radius contracts down to a ZAMS radius 


5 


Main Sequence 




Star has contracted enough for Tc to reach ~ 10^ K 
Hydrogen ignites and star stabilizes onto the main sequence. 



3 STAR FORMATION AND FEEDBACK IN 
THE CLUSTER ENVIRONMENT 

We investigate what difference protostellar modeling makes 
to the overall evolution of the cluster. The most important 
consequence of the improved hybrid characteristics raytrac- 
ing code employed by Peters et al. (2010a) is that it allows 
for the realistic simulation of HIT regions, with ionization, 
heating, and shadowing effects built in. One of the most im- 
portant consequences of accurate pre-main-sequence model- 
ing is that it tempers the ionization and heating in the early 
stages of star formation. 

To study this effect, we look at two variables: mean ion- 
ization in our simulation box, and mean gas temperature. 
Figure [4] shows these two measurements as functions of time 
in our simulation. Time is measured in units of global freefall 



time, or t^^ ^ 590, 000 years. With the model of lPeters et al.l 
f2010a), sink particles follow a ZAMS model for the stel- 
lar radius and intrinsic luminosity, which means that they 
are hotter and more compact than true pre-main-sequence 
stars. This causes them to release more ionizing photons, 
compared to the protostellar case. The onset of ionization 
in this case leads the evolving protostar case by about 0.03 
freefall times, or about 17.7 kyr. The onset of star forma- 
tion in both simulations occurs at around 1 freefall time. 
For both cases, after 1.1 freefall times, the largest star in ei- 
ther simulation is at ^ 2OM0 and dominates the UV output 
of the cluster, resulting in comparable mean ionization for 
both cases. 

When we consider mean temperature instead of mean 
ionization, the leading effect by the ZAMS-based model is 
still there, only less pronounced. Major heating of the gas 
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Figure 4. A comparison of the mean ionization fraction and mean 
temperature in cluster simulations with different prestellar mod- 
els. In each case, the mean is calculated by finding the volume- 
weighted average. Values are only meaningful in a relative sense, 
as the simulation volume is large (side length ~ 3.8pc) and the 
most active region is the inner cubic parsec. 



in this case leads the evolving protostar model case by close 
to 0.01 freefall times, or about 5.9 kyr. The first star to 
form in a cluster tends to grow to be among the largest 
stars in the cluster and dominate the heating and ionization. 
This suggests that accurate protostellar modeling is most 
important in th e early stages of a cluster simulation, and for 
low-mass stars. I Off ner et al.l (2009) showed that radiation 
even from low-mass stars has a significant effect on the gas 
heating and formation of brown dwarfs. 

To get a visual sense of the gas dynamics and configu- 
ration of the cluster, we visualize the gas density by taking 
slices through our simulation box. Zoomed-in views of the 
cluster are shown in Figure O The simulation box is actu- 
ally about 3.8 pc across. Here we show the central region, at 
about 0.5 pc across. The upper row in the figure shows the 
simulation results with the evolving protostar model, while 
the lower row shows the A-ZAMS results. In each row, the 
panels show: gas density in a horizontal slice through mid- 
plane of the simulation box (left), gas density in a vertical 
slice showing the cluster edge-on (centre), and gas temper- 
ature in the same vertical slice (right). Gas temperature is 
discussed in Section 13.41 The two density panels also show 
velocity vectors for the high- velocity gas. The fastest-moving 
gas travels at close to 30 km/s. Gas densities range from 
10"^^ to 10"^^ g/cm^. The hollowed-out HII regions, where 
the gas is largely ionized, expand outwards above and below 
the disc as a kind of fountain before falling back onto the 
disc. 



Sink particles indicating the locations of stars are 
marked with black-rimmed gray points. The side view shows 
the stars to be confined to the disc while the top-down shows 
the stars packed in a tight cluster. The separation between 
stars nowhere exceeds 0.1 pc. During the course of the sim- 
ulation, stars are seen to be dynamically interacting, ex- 
changing angular momentum, forming and breaking apart 
binaries. 

These snapshots of the simulation are taken at around 
1.21 freefall times in each case, near the end of the simu- 
lation. At this stage, about 714 kyr have elapsed since the 
beginning of the simulation, with the onset of star formation 
having occurred at around 600 kyr. At this stage, both model 
results look similar in many ways: the stars are in a densely- 
packed cluster, and each cluster has produced an expanding 
HII region. The HII regions in each figure are approximately 
the same size, although amorphous and variable. They do 
not seem to be affected by our choice of prestellar model. 
This is because of how each cluster has become dominated 
by massive stars already evolved onto the main sequence, 
and the differences between prestellar models has vanished. 
The stars are all releasing copious amounts of ionizing radi- 
ation, driving the evolution of these HII regions. 



3.1 Binaries 

IZinnecker Yorkd (|2007l l state that massive stars occur 
more frequently in binaries relative to low-mass stars. Lack- 
ing turbulence and magnetic fields, our molecular gas clumps 
do not represent the true initial conditions for cluster for- 
mation, but the stars in our simulations to form binaries. 
There is no reason to suspect that the choice of protostellar 
model has any effect on binary formation or binary mass ra- 
tios. Lacking ensemble averages, we cannot make any special 
claims, but report that of the 5 stars formed in each of our 
simulation, 4 stars end up in binaries. Dynamical encoun- 
ters between stars cause binaries to form, break apart, and 
reform. The final mass ratios of the two pairs each simula- 
tion were 3.74 and 1.61 with the Offner model. The Peters 
model simulation saw mass ratios of 3.38 and 1.92. The final 
masses of the stars formed in each simulation are reported 
in Tabled 



3.2 Accretion Histories 

We now compare the simulations wit h a focus on the ac - 
cretion histories of the sink particles. IPeters et al ] (|2010al l 
have shown that the gas surrounding the centre of the clus- 
ter would fragment and result in a highly variable accre- 
tion rate. We see this in Figures [6] and [71 where we show in 
the two panels the accretion histories of every sink particle 
formed in our simulation along with their accretion rates. 

In the lower panel we see the accretion rate of each star, 
and for most of the stars in our simulation, the accretion rate 
remains between lO"'^ and Mq /yi. 

The upper panel in Figure [6] shows the growing masses 
of each of the stars in our protostellar model simulation. 
Star formation does not really commence until after the first 
dynamical time (freefall time) — about 0.59 Myr for our sim- 
ulation setup. There seems to be a burst of star formation 
after t ~ l.lOt^. Interestingly, the most massive star is not 
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Densily (g/cm'^3) 

l.Oe-15 



l.Oe-19 



4.6e-2 



Density (g/cm-^S) 

l.Oe-15 



l.Oe-19 



4.6e-2 



2.2e-22 



Temperature <K) 

20000, 



6766.7 



88.669 



30.000 




774.60 



30.000 



Figure 5. The two rows show the results fr om the two prestella r models tested, with the A-ZAMS model of IPeters et al.l (|2010al ) in the 
top row and the evolving protostar model of lOffner et alJ (|2009l ) in the bottom row. Each is shown near the end of the simulation, after 



about 1.21 freefall times (714 kyr). In each row, the panels show, from left to right: the gas density in a horizontal slice through the 
midplane, the gas density in a vertical slice in the centre of the simulation box, and the gas temperature in the same slice. Scale bars 
indicate the physical sizes and the speeds represented by vectors in the gas density panels. The scale for these vectors is the same for 
both side views and top-down views. Stars are indicated by black circles. 



the first star in our simulation, but it is overtaken in mass by 
the second star, which reaches a final mass of about 43.5M0. 
The others reach final masses of approximately 36.0, 28.5, 
22.3, and 7.6 Mq. The average mass of these 5 stars is 27.6 
Mq. We were able to run the evolving protostars simula- 
tion longer than the A-ZAMS case. During this extra time, 

3 additional stars formed and accreted about IMq of ma- 
terial each, but we do not use this additional data in our 
comparison with the A-ZAMS case. 

In the A-ZAMS case, shown in Figure[7l the final masses 
of the stars are 47.3, 33.3, 28.4, 17.3, and 14.0 Mq. The av- 
erage mass of these 5 stars is 28.1 Mq. It is difficult to draw 
firm conclusions about the impact of a protostellar model 
on the population dynamics of a cluster. We would need to 
complete longer simulations under more realistic initial con- 
ditions (including turbulence). The evolving protostars run 
experienced a second wave of star formation, but when we 
restrict ourselves to comparing only the first 1.21 in each 
case, we find that they have almost the same average mass. 

Interestingly, though, the evolving protostars case had 

4 stars with masses greater than IOM0. These were all more 



Table 3. Final stellar masses after 1.21 freefall times (714 kyr) 
in each cluster simulation, in units of solar masses, comparing the 
different prestellar model results. 



Evolving Protostars 


Augmented ZAMS 


43.5 


47.3 


36.0 


33.3 


28.5 


28.4 


22.3 


17.3 


7.6 


14.0 



closely packed (smaller variance), than the four most mas- 
sive stars in the A-ZAMS case. We propose that the reduced 
initial heating and ionization from the self-consistently 
evolved pre-main-sequence stars results in a more equitable 
partition of mass between the massive stars. The most mas- 
sive star in this simulation outranks the second largest by 
about 7.5M0. By comparison, the leading star in the A- 
ZAMS simulation exceeded the next most massive star by 
14M0, or nearly double. Further simulations with different 
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Time (t^) 





Time (t^) 
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Figure 6. Accretion histories of stars formed in the cluster sim- 
ulation of the evolving protostars setup. The upper panel shows 
the mass of each particle as a function of time. The lower panel 
shows the accretion rate in units of Mq yT~^ as a function of 
time. The dynamical time is about 0.59 Myr. 



Figure 7. Accretion histories of stars formed in the cluster sim- 
ulation of the A-ZAMS setup. The upper panel shows the mass 
of each particle as a function of time. The lower panel shows the 
accretion rate in units of Mq yT~^ as a function fo time. The 
dynamical time is about 0.59 Myr. 



initial conditions are required to confirm whether this choice 
of prestellar model will always have such an impact. 

These are the results of only a single simulation in each 
case, so it is difficult to say that this difference in massive 
spectrum is highly significant, especially given that the av- 
erage mass of each cluster is similar and our simulations did 
not contain turbulence. By other measures, such as the av- 
erage ionization, mean temperature, or HII region morphol- 
ogy, the two prestellar models converged and gave similar 
results. The mass spectrum shows a similar average mass of 
28Mq , but with the evolving protostar model having both 
a smaller peak mass and smaller difference in mass between 
the two top stars relative to the A-ZAMS model. 

We were able to run the protostellar simulation a little 
longer than the A-ZAMS simulation; there occurred a sec- 
ond burst of star formation that only appeared very late in 
the simulation. These stars grew to be 1.9, 1.3 and I.OMq. 
The A-ZAMS run may have formed more stars if run for 
longer. We have run each setup for approximately two weeks 
on 64 processors, or approximately 21,500 CPU- hours. The 
protostellar simulation progressed further than the ZAMS 
simulation. In either case, memory or eventual code stabil- 
ity limited the length of the runs. 

3.3 Mass-radius relation 

The mass-radius relation for a star is a means of compar- 
ing different protostellar models. It is also a way of seeing 



the evolution of the stars in our simulation. As stars accrete 
mass or undergo nuclear-structural changes in their interi- 
ors, the radius reacts either by expanding or contracting. We 
see the evolution of the stars in our simulation represented 
in Figure [8] In this figure we compare the radii of stars from 
our different cluster simulations. 

We show the accretion radius for a single star in the 
A-ZAMS run by the gray line in Figure [8] The red lines in 
this figure are the stars following the protostellar evolution 
model that we have described in Section 12.11 These stars 
have their radius continuously evolved according to their 
burning state and the accretion of new material. The ra- 
dius, therefore, does not fluctuate with unrealistic rapidity. 
Because the model has a self-consistent description of the 
radius, we use the same quantity to describe the stellar ra- 
dius and the accretion radius, rather than computing each 
by different means. Protostars have radii an order of mag- 
nitude larger than a zero-age main-sequence star of equal 
mass. Hence, their effective temperatures and flux of ioniz- 
ing photons are going to be much less (for a IMq star, 3000K 
vs. 5000K in effective temperature, 10^^ s~^ vs 10^^ s~^ in 
ionizing photons). Star in simulations without protostellar 
modeling may excessively heat or ionize the gas during the 
early phases of star formation. 

The evolution of the stars in each simulation is also 
revealed by the mass-luminosity relation, shown in Figure 
El Black and grey lines belong to a representative star in 
the the A-ZAMS model simulation, dark red and blue to a 
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ZAMS 
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in Peters et al. (2010a) 







Mass {Mq) 

Figure 8. The mass-radius relation for the stars in both cluster 
simulations. The black dashed line marks the stellar radius track 
of a sink particles in the A-ZAMS simulation. The radius is based 
on tabulated values of luminosities and temperatures for ZAMS 
stars. The gray line indicates the separately-calculated accretion 
radius for a single star. Red lines mark the tracks of sink particles 
following the evolving protostar model. This protostellar radius 
is used as both the stellar radius and the accretion radius. 

representative star in the protostellar simulation. Accretion 
luminosity, calculated as Lace = GMM/Racc, is especially 
sensitive to the accretion rate M and the accretion radius 
Race- Stars show accretion luminosities that are up to an 
order of magnitude larger than the stars in the evolving 
protostar model simulation on account of the difference in 
stellar radius. Only for stars larger than about 20Mq do the 
differences between the two models disappear. The black 
jagged line indicates the main sequence luminosities from 
a precomputed table, which tends to underestimate stellar 
luminosities for protostars less than about 3Mq . Protostellar 
luminosity of one star is given by the dark red line and 
protostellar accretion luminosity by the blue line. Much of 
the rapid fluctuation in the accretion luminosities of both 
simulations stems from the highly variable mass accretion 
rate (see Figures [6] and [T)) . 

3.4 Ionization and Temperature 

Protostars have large radii about an order of magnitude 
larger than equivalent- mass main-sequence stars. They may 
be just as luminous, and they certainly have high accretion 
luminosities, but it is the effective temperatures of their sur- 
faces that determine how great the flux of ionizing photons 
will be, if the star emits any at all. The single greatest differ- 
ence we see when simulating the evolution of a star cluster 
with self-consistent protostellar modeling is that when the 
first stars begin to form after about a dynamical time, the 
average gas temperature and average ionization of the gas is 
considerably less in the simulation involving our protostellar 
model (Figure 

Figure [TOl shows the mass- weighted spectrum of the ion- 
ization fraction in both cases, with the evolving protostar 
model on the left and the A-ZAMS model on the right. 
Values for ionization fraction range from 10"^*^ to approx- 
imately 1 (completely ionized). Ionization fractions <^ 1 
should not be taken too seriously, as our model includes 
only stellar ionizing radiation from the stars in our cluster. 
The figure shows the spectrum for all the gas involved in 




Mass {Mq) 

Figure 9. The mass-luminosity relation for a representative star 
in each of the two cluster simulations. For each star, its accre- 
tion luminosity and intrinsic stellar luminosity are plotted. The 
red and dashed blue lines show the instrinsic stellar luminosity 
and the accretion luminosity, respectively, of a star in the evolv- 
ing protostars simulation. The black dashed line, with its stepped 
appearance, represents the intrinsic luminosity of a ZAMS star, 
retrieved from a table of ZAMS values. Finally, the grey line shows 
the accretion luminosity of a star in our in the A-ZAMS simula- 
tion. The luminosity is calculated as in Pe ters et al. (2010a) by 
an interpolation of the radius to models bv lHosokawa Sz, Omukail 
(■200ftV 



the simulation — approximately 1000 Mq in total. The thick 
yellow line indicates the mass- weighted average value for the 
ionization fraction with the value printed beside the line. In- 
dividual snapshots in time are: t = 1.05, 1.10, 1.15, 1.20 t^jp. 

It is important to show how the averages change over 
time in Figure |4] because of how the mean tends to fluctuate 
yet the two models have similar values for all but the earliest 
phases of star formation. The early phase is shown in the 
first row, at t = 1.05t^. Here there is a significant difference 
in the mean ionization fractions of the two models. The low- 
mass ZAMS stars are hotter and have smaller radii. There is 
greater early ionization seen in this case. At later stages, the 
distributions appear more similar as the conspicuous effects 
of the model disappear. 

The temperature structure of the gas surrounding the 
cluster is shown in right panels of Figure [5] These reveal 
some interesting features. The gas in the vicinity of the 
cluster is approximately 100 K, heating to this tempera- 
ture by the nonionizing radiation coming from the clus- 
ter. We also see pockets of very hot (10^ K), ionized gas 
in the expanding HII regions. When we study the evolu- 
tion of these regions in time, we see that these pockets of 
high-temperature gas are very transient: forming, expand- 
ing, breaking apart, and cooling very rapidly. They are due 
to photoionization and photoion ization heating caus ed by 
the massive stars in the cluster. IPeters et al.l (|2010al ) have 
attributed this flickering to the chaotic gas motions in the 
cluster. Gas moves toward the interior of the cluster through 
the disc and interacts with the ionizing radiation giving 
rise to many unstable morphologies that expand outwards 
above and below the disc. This has the appearance of flicker- 
ing on relatively short timescales: less than 560 years — the 
temporal resolution of our simula t ions. S ynthetic observa- 
tio ns of the original Peters et al.l (jloiOal) results analyzed 
bv lGalvan-Madrid et al.1 (|201lh have revealed this flickering 



© 2011 RAS, MNRAS 000, [THE] 



10 



M. Klassen, R.E. Pudntz, & T. Peters 



Offner et al. (2009) Model 



Peters et al. (2010a) Model 




-10 -8 -6 -4 -2 
Log Ionization Fraction 



10 



-8 -6 -4 -2 
Log Ionization Fraction 



Figure 10. The evolving mass -weighted ioniz ation fraction spectrum. Compared a re cluster simulation s with stars running on the 
evolving protostar model of O ffner et al.l (|2009i ') on the left and the A-ZAMS model of I Peters et al.l (l2010a| ) on the right. A distribution 
of the total mass in the simulation box (about 1000 Mq) is shown for t = 1.05, 1.10, 1.15, 1.20 1^. One freefall time is approximately 0.59 
Myr. The yellow line indicates the mass-weighted average ionization fraction, the numerical value of which is printed to the left of the 
line. 



visible in radio-continuum emission and have demonstrated 
that it is in agreement with available observations. 



4 DISCUSSION 

Considering that the radii and luminosities of true pro- 
tostars are vastly different from their ZAMS counterparts 
(i.e. Figures [8] and (9]) , it may seem surprising that our two 
simulations actually look so alike. For instance, the two sim- 
ulations form an equal number of stars, their average mass is 
approximately the same, and morphology of the clump with 
its outflows and HII regions appear qualitatively similar in 



both cases. It is important to note what we are com par- 
ing. T he model that we are comparing against (Peters et al.l 
l2010al ) treated the stellar radius and the accretion radius 
separately, meaning that gas heating has two components: 
one due stellar radiation, and one due to the accretion lumi- 
nosity. The mass- luminosity relation of Figure [9] shows that 
the ZAMS stellar luminosity underestimates the true proto- 
stellar luminosity for pre-main-sequence stars. It also shows 



that the accretion luminosity, calcula ted as injPeters et al 
(2010a) by an interpolation to the iHosokawa Omukai 



(2009) models, overestimates the true protostellar accretion 
luminosity. So there are two competing differences and these 
effects will partially cancel each other out. The result is that 
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our evolving protosta rs simulation looks simi lar in many 
ways to the results of IPeters et alj (|201Qal V If IPeters et alj 
had not boostrapped the separate treatment of ac- 
cretion radius on to the ZAMS model, there may have been 
a gross overestimation of the accretion luminosity — which 
dominates the total luminosity of a star during its early life- 
time. The errors resulting from this overestimation could be 
substantial. 

In our radiative feedback technique, we treat ionization 
separately from heating, and ionization depends solely on 
the effective temperatures of our stars. Since protostars have 
cooler surface temperatures than ZAMS stars of equal mass, 
there is much less early ionization. Since it is the ZAMS stars 
of high mass that dominate the radiation output of a clus- 
ter, the differences between our model and the ZAMS model 
dissapear after the early stages of stellar evolution (Figure 
nop . A side- by-side comparison the ionizing flux from stars 
with different stellar models will be included in a forthcom- 
ing paper. 

Our simulations have a number of limitations that 
should be noted. They neglect the effects of radiation pres- 
sure. On large scales, radiation pressure from stell ar clus- 
ters could drive galactic winds (jMurrav et ahllioilh . How- 
ever, within our low-density lOOOM© cluster, radiation pres- 
sure below the Eddington limit should not be dynamically 
signif icant |Yorke Sonnhalterll2002l : lKrumholz Matznerl 
I2OO9I I. After the first absorption/reemission event, the radi- 
ation will have been converted to infrared radiation to which 
the molecular cloud is largely transparent. The first absorp- 
tion event is unlikely to impart a significant amount of mo- 
mentum. 

In our simulations, we have treated gas that was initially 
cold and in solid body rotation, but without any turbulence. 
Cluster-forming clumps in molecular clouds are observed to 
have supersonic turbulence, and a more realistic set of initial 
conditions would include turbulence. However, this might 
have obscured the effects of our protostellar model that we 
were seeking to measure. We are currently preparing to run 
simulations that include realistic turbulent initial conditions 
as well magnetic fields , which were also le ft out of this sim- 
ulation (see, however, IPeters et al.l (|201lh for the effects of 
magnetic fields on our non-turbulent initial conditions). 

The protostellar model we have added to our simula- 
tions improves on previous work by adjusting the ionizing 
luminosity so that it matches the stellar surface effective 
temperatures for accreting protostars, which initially have 
radii larger than equal-mass stars on the main sequence. 
We note, however, that a full-spectrum treatment of the ra- 
diation still faces technical and computational limits that 
make the problem extremely challenging. As a compromise, 
we break the radiation into its ionizing and nonionizing com- 
ponents. 



5 CONCLUSIONS 

Stars begin to affect their birth environments as soon as they 
are born through radiative feedback. We have considered 
the impact that pre-main-sequence modeling can have on 
a star cluster by comparing two different prestellar models 
already described in the literat u re. We did this by repeating 
the simulation of IPeters et al.l (|2010al V We then upgraded 



the FLASH code to include a protostellar evolutio n module 
based on the one described in the appendices of .Offner et al.l 
([2OO9I ). 

Each model works by equipping the stars in the sim- 
ulation ("sink particles") with a stellar radius and lumi- 
nosity. The greatest difference between th e two models was 
self-consistency. The Peters et al.| (|2010al l model calculated 
approximate stellar parameters on-the-fiy, while the evolv- 
ing protostar model evolved the stellar parameters self- 
consistently through the simulation as the stars grew and 
accreted mass. 

In terms of the overall gas structure, HII regions, tem- 
perature structure, mean ionization fraction, or stellar bina- 
rity, the two models produced qualitatively the same results. 
This is because a cluster comes to be dominated by its most 
massive stars, which are evolved, main-sequence, highly lu- 
minous stars, regardless of the choice of stellar model. These 
one or two massive stars control the overall dynamics. 

The differences exist in the early phase of star forma- 
tion. Major ioniz ation of the gas in t he evolving protostar 
model lagged the lPeters et al.l (|2010al l model by about 3% 
of a freefall time, or about 17.7 kyr. Major heating of the gas 
lagged by about 1% of a freefall time, or about 5.9 kyr. The 
difference in heating and ionization was due to the fact in 
Peters et al. 1 (|2010al \ the stellar radius was underestimated 
(a Z AMS-equivalent value was taken) , when protostars have 
radii an order of magnitude larger than a zero-age main- 
sequence star of equal mass. The correspondingly higher sur- 
face temperatures resulted in excess heating and ionization 
in this model. When both models had stars converging onto 
the main sequence, the differences between the two models 
diminished. 

It is possible that these initial differences could have had 
a lasting effect on the stellar population. The most massive 
star at the end of each simulation was 43.5M(:) in the evolv - 
ing protostar model, and 47. SM© in the Pe ters et al.l (|2010ah 
case — a difference of 8%. The differences in mass between 
the most massive star and the next largest star_was_7^5Mo 
i n the e volving protostar case and 14M0 in the lPeters et al.l 
(|2010al ) case. It would require further simulations, varying 
the initial conditions, to confirm that this is always the case. 

The cluster of stars is embedded in a rotating disc of 
gas approximately 0.2 pc in size. The expanding HII regions 
above and below the disc are rapidly changing in shape and 
size on timescales shorter than 570 years. The physical size 
of these HII regions in our simulation is at most about 0.2 
pc. This flickering is observed regardless of the prestellar 
model used. 

Future simulations will have initial conditions includ- 
ing turbulence to model molecular clouds as realistically as 
possible. The stars will no longer be forming within a global 
disc, but rather along sheets and filaments in diverse parts 
of the cloud. With star formation thus spread out more in 
space and time, we expect the influence of individual young 
stellar objects on their environments to be more significant 
than when all stars form in a central cluster. It will be im- 
portant to have the radiative feedback accurately modeled 
in these cases. 



© 2011 RAS, MNRAS OOO.fTIO 



12 M. Klassen, R.E. Pudntz, & T. Peters 



ACKNOWLEDGMENTS 

We thank Mark Krumholz for clarifying some of the more 
technical aspects of the protostellar evolution method em- 
ployed in the ORION code and described in tOffner et al.l 
(|2009l l. This was very helpful in developing the FLASH mod- 
ule following the same method. We thank Takashi Hosokawa 
for sharing data from Hosoka wa Omukail (|2009') with us. 
We also thank our anonymous referee for a very careful re- 
view of our paper that helped to significantly clarify the 
presentation of our results. M.K. acknowledges financial 
support from the Ontario Graduate Scholarship Program. 
R.E. P. is supported by a Discovery Grant from the Nat- 
ural Sciences and Engineering Research Council (NSERC) 
of Canada. T.P. acknowledges financial support as a Fellow 
of the Baden- Wiirttemberg Stiftung funded by their pro- 
gram International Collaboration II (grant P-LS-SPII/18) 
and through SNF grant 200020_137896. The FLASH code 
was in part developed by the DOE-supported Alliances Cen- 
ter for Astrophysical Thermonuclear Flashes (ASCI) at the 
University of Chicago. This work was made possible by 
the facilities of the Shared Hierarchical Academic Research 
Computing Network (SHARCNET: www.sharcnet.ca) and 
Compute/Calcul Canada. 



REFERENCES 

Banerjee R., Pudritz R. E., Anderson D. W., 2006, MN- 

RAS, 373, 1091 
Banerjee R., Vazquez- Semadeni E., Hennebelle P., Klessen 

R. S., 2009, MNRAS, 398, 1082 
Bate M. R., 2009, MNRAS, 392, 1363 
Beuther H., Churchweh E. B., McKee C. F., Tan J. C, 

2007, Protostars and Planets V, 165 

Blitz L., 1993, in Protostars and Planets HI, E. H. Levy & 
J. I. Lunine, ed., pp. 125-161 

Clarke C. J., Bonneh I. A., Hillenbrand L. A., 2000, Pro- 
tostars and Planets IV, 151 

Evans, II N. J., 1999, ARA&A, 37, 311 

Federrath C, Banerjee R., Clark P. C, Klessen R. S., 2010, 
ApJ, 713, 269 

Franco- Hernandez R., Rodriguez L. F., 2004, ApJL, 604, 
L105 

Fryxeh B. et al., 2000, ApJS, 131, 273 

Galvan-Madrid R., Peters T., Keto E. R., Mac Low M.-M., 

Banerjee R., Klessen R. S., 2011, MNRAS, 1081 
Galvan-Madrid R., Rodriguez L. F., Ho P. T. P., Keto E., 

2008, ApJL, 674, L33 

Gomez L., Rodriguez L. F., Loinard L., Lizano S., Allen 
C, Poveda A., Menten K. M., 2008, ApJ, 685, 333 

Hoare M. G., Kurtz S. E., Lizano S., Keto E., Hofner P., 
2007, Protostars and Planets V, 181 

Hosokawa T., Omukai K., 2009, ApJ, 691, 823 

Iliev I. T. et al., 2006, MNRAS, 371, 1057 

Keto E., 2002, ApJ, 580, 980 

Keto E., 2003, ApJ, 599, 1196 

Keto E., 2007, ApJ, 666, 976 

Krumholz M. R., Cunningham A. J., Klein R. I., McKee 

C. F., 2010, ApJ, 713, 1120 
Krumholz M. R., Klein R. I., McKee C. F., 2007, ApJ, 656, 

959 



Krumholz M. R., Klein R. I., McKee C. F., 2011, ApJ, 740, 
74 

Krumholz M. R., Matzner C. D., 2009, ApJ, 703, 1352 
Matzner C. D., 2002, ApJ, 566, 302 
McKee C. F., Tan J. C, 2003, ApJ, 585, 850 
Mezger P. G., Henderson A. P., 1967, ApJ, 147, 471 
Murray N., Menard B., Thompson T. A., 2011, ApJ, 735, 
66 

Nakano T., Hasegawa T., Morino J.-L, Yamashita T., 2000, 
ApJ, 534, 976 

Nakano T., Hasegawa T., Norman C, 1995, Ap&SS, 224, 
523 

Offner S. S. R., Klein R. I., McKee C. F., Krumholz M. R., 

2009, ApJ, 703, 131 
Palla F., Stahler S. W., 1991, ApJ, 375, 288 
Palla F., Stahler S. W., 1992, ApJ, 392, 667 
Paxton B., 2004, PASP, 116, 699 

Peters T., Banerjee R., Klessen R. S., Mac Low M.-M., 

2011, ApJ, 729, 72 
Peters T., Banerjee R., Klessen R. S., Mac Low M.-M., 

Galvan-Madrid R., Keto E. R., 2010a, ApJ, 711, 1017 
Peters T., Klessen R. S., Mac Low M.-M., Banerjee R., 

2010b, ApJ, 725, 134 
Peters T., Mac Low M.-M., Banerjee R., Klessen R. S., 

Dullemond C. P., 2010c, ApJ, 719, 831 
Pollack J. B., Hollenbach D., Beckwith S., Simonelh D. P., 

Roush T., Fong W., 1994, ApJ, 421, 615 
Rijkhorst E.-J., Plewa T., Dubey A., Mellema G., 2006, 

A&A, 452, 907 
Rodriguez L. F., Gomez Y., Tafoya D., 2007, ApJ, 663, 

1083 

Spitzer L., 1978, Physical processes in the interstellar 

medium. New York Wiley-Interscience, 1978. 333 p. 
Tan J. C, McKee C. F., 2004, ApJ, 603, 383 
Testi L., Sargent A. I., Olmi L., Onello J. S., 2000, ApJL, 
540, L53 

Tout C. A., Pols O. R., Eggleton P. P., Han Z., 1996, MN- 
RAS, 281, 257 

Williams J. P., Blitz L., McKee C. F., 2000, Protostars and 

Planets IV, 97 
Wood D. O. S., Churchweh E., 1989, ApJS, 69, 831 
Yorke H. W., Sonnhalter C, 2002, ApJ, 569, 846 
Zinnecker H., Yorke H. W., 2007, ARA&A, 45, 481 



© 2011 RAS, MNRAS 000,[THl2] 



